function output = runNeuralNet(net, in )
    f = net.f;  %function handle
    W = net.W;  %list with matrices for each layer k: W{k}
    L = length(W); %nr of layers
    
    input{1} = [in; 1];
    for k=1:L %heen
        activation{k} = W{k}*input{k};
        output{k} = f(activation{k});
        input{k+1} = [output{k}; 1]; %add bias
    end
    output = output{L};
end

